R code and downstream analysis objects for the scRNA-seq atlas of normal and tumorigenic human breast tissue

Breast cancer is a common and highly heterogeneous disease. Understanding cellular diversity in the mammary gland and its surrounding micro-environment across different states can provide insight into cancer development in the human breast. Recently, we published a large-scale single-cell RNA expression atlas of the human breast spanning normal, preneoplastic and tumorigenic states. Single-cell expression profiles of nearly 430,000 cells were obtained from 69 distinct surgical tissue specimens from 55 patients. This article extends the study by providing quality filtering thresholds, downstream processed R data objects, complete cell annotation and R code to reproduce all the analyses. Data quality assessment measures are presented and details are provided for all the bioinformatic analyses that produced results described in the study.

Next, breast tissue from BRCA1 mutation carriers was obtained for investigating cellular changes in precancerous state. Overall, the differences of stromal and immune subsets between normal and BRCA1+/− preneoplastic tissue were not significant, nor was the proportions of different cell clusters. However, extensive changes in the tissue micro-environment were observed between the preneoplastic and the neoplastic states in BRCA1 mutation carriers 5 .
Finally, ER+, HER2+ and triple negative breast cancer (TNBC) tumors were obtained from treatment-naive patients for exploring the degree of heterogeneity within the cancer cell compartment and its micro-environment across different tumor subtypes. Extensive inter-patient heterogeneity was revealed by single cell integration analyses across all cancer subtypes. Within the tumor populations, a discrete cluster of cycling MKI67+ tumor cells were observed for all three major breast cancer subtypes. Within the tumor micro-environment, different immune landscapes were observed in different cancer subtypes. Both TNBC and HER2 featured a proliferative CD8+ T-cell cluster, whereas ER+ tumors primarily comprised cycling TAMs. In addition, matched pairs of ER+ tumors and involved lymph nodes were profiled for examining the relationship between primary breast tumors and malignant cells that seed lymph nodes. Clonal selection and expansion were observed in some patients, whereas mass migration of cells from the primary tumor to the LN was observed in some other patients 5 .
The ScBrAtlas provides a valuable resource for understanding cellular diversity and cancer genesis in human breast. The examination and exploration of the single cell data presented in this study required large-scale bioinformatics analyses for multiple groupings of the original data. While genewise read counts were previously made

Male ER+ tumor (2)
In this report we describe the bioinformatics analysis used in the ScBrAtlas in greater detail. We provide a complete description of the quality control filters used to select 341,874 cells for downstream analyses. The technical quality of both the 10X single-cell transcriptomic data sets and the bulk RNA-seq reference data set is assessed to demonstrate the reliability of the data. We provide downstream R data objects corresponding to each data integration and cell clustering presented in the ScBrAtlas, together with R code to reproduce the data objects. Crucially, the data objects provided here include cell barcodes by which each individual cell can be tracked through all the analyses. We also provide detailed information allowing the copy number variation analyses to be mapped back to individual samples and cell clustered, thus providing a way to distinguish putative malignant cancer cells from normal epithelial cells in the cancer tumors. All the resources and the detailed information can be easily accessed and utilized by researchers for further exploration and clinical validation, which may lead to discoveries of novel approaches for personalized breast cancer treatment in the future.

Methods
Human Samples. This article is a companion to the ScBrAtlas study and ethics approval regarding human samples is as stated in that article 5 . Human breast tissues were obtained from consenting patients through the Royal Melbourne Hospital Tissue Bank, the Victorian Cancer Biobank and kConFab with relevant institutional review board approval. Human Ethics approval was obtained from the Walter and Eliza Hall Institute Human Research Ethics Committee.

Read alignment and count quantification.
Single-cell RNA-seq expression profiles of 69 samples from 55 patients were generated using the 10x Genomics Chromium platform and an Illumina NextSeq 500 sequencer (Fig. 1a, Supplementary Table 1). Genewise read counts were produced for all samples using Cell Ranger v3.0.2 (https://support.10xgenomics.com). Specifically, the Illumina base call files (BCLs) were demultiplexed into FASTQ files by "cellranger mkfastq" then genewise read counts were obtained using "cellranger count" with the Cell Ranger human GRCh38 reference v3.0.0. Default settings were used for all parameters apart from file locations and memory usage. For each biological sample, results were taken from the Cell Ranger output directory "outs/filtered_feature_bc_matrix". The output for each sample consists of three files: the count matrix in matrix market mtx.gz format, the cellular barcodes (barcodes.tsv.gz) and the gene identifiers (features.tsv.gz). All samples share the same gene identifiers but the count matrix and barcode files are sample-specific (Supplementary Table 1). The files provide results for a total of 421,761 cells across the 69 samples (Supplementary Table 2). All cells included in the Cell Ranger output have least 500 reads successfully assigned to genes. The Cell Ranger output files were deposited as GEO series GSE161529 7 and were used for downstream bioinformatics analyses.
Quality control and cell filtering. All samples were first checked using marker genes to confirm the presence of epithelial cells and the absence of significant stromal cell contamination. Then an average of 15% of cells were filtered from each sample based on (i) the total number of mapped reads for that cell (library size), (ii) the number of genes detected and (iii) the proportion of reads mapping to the mitochondria 8 (Fig. 2). A lower bound of 500 was generally applied to the number of detected genes for each cell, although this was reduced to 400 or 300 for a small number of samples with low read coverage. Upper bounds were applied to the number of detected genes and to library size to remove potential doublets. Thresholds were chosen for each sample by plotting genes detected versus library size and choosing thresholds that excluded outliers but kept the main body of cells.
For most samples, an upper limit of 0.2 was placed on the mitochondrial proportion with 80-100% of cells below this threshold. For one high quality sample with very low mitochondrial proportions the threshold was tightened to 0.1. For a few samples with higher mitochondrial proportions, the mitochondrial threshold was loosened to 0.25, 0.3 or 0.4 in order to keep the majority of the cells. These samples were considered of lesser quality but not so poor as to be excluded from the analysis.
The threshold values used for these QC metrics are shown in Supplementary Table 2 and are also supplied in machine-readable form as part of the data submission 9 . A total of 341,874 cells remained after quality filtering for downstream analysis.
Single-cell RNA-seq integration analysis. The samples included breast tissues from normal healthy donors, BRCA1 mutation carriers and patients diagnosed with different types of breast cancer (triple negative, ER+ and HER2+). Matching pairs of tumor and lymph node (LN) samples, as well as tumor samples from male patients, were also included. The single-cell analysis strategy involved grouping together comparable samples, integrating the profiles, then clustering cells into putative cell types. A total of 16 different sample-groups were integrated (Fig. 1b). Some samples were involved in more than one integration, for example the pre-neoplastic samples with BRCA1 mutations were integrated first with the normal samples and later with the BRCA1 triple negative (TN) tumor samples. For some sample-groups analyses, subsets of cells were extracted, re-integrated and re-clustered. The total number of cell cluster analyses is shown in Table 1.
Samples were integrated using the Seurat anchor-based integration method 10 . To perform dimensionality reduction, the first 30 principal component were computed and used for the cell clustering and t-distributed stochastic neighbor embedding (t-SNE) visualization 11 . The default Louvain clustering algorithm 12 was used for cell cluster identification. Cluster resolutions were typically set to values around 0.1, lower than the Seurat default, in order to ensure conservative and reproducible clusters. The only exception was Figure EV4A where much higher resolutions were used in order to distinguish subsets of T cells 5 . (2022) 9:96 | https://doi.org/10.1038/s41597-022-01236-2 www.nature.com/scientificdata www.nature.com/scientificdata/ We provide here the Seurat data objects containing each of the cluster analyses as R data files ( Table 2). The R data objects contain cell cluster details for each cell. The R code by which each R object was constructed is also provided ( Table 2). Differential expression and pathway analysis. Differential expression analyses were performed to detect marker genes for different cell clusters. In order to account for the biological variation between different patients, a pseudo-bulk approach was used in most cases where read counts from all cells under the same cluster-sample combination were summed together to form pseudo-bulk samples. The edgeR's quasi-likelihood pipeline was used for pseudo-bulk differential expression analysis, where the baseline differences between patients were incorporated into the linear model 13 . The Seurat's FindMarkers function was applied where pseudo-bulk samples were not satisfactory due to low cell numbers or imbalanced cluster-sample combination. KEGG pathway analyses were performed using the kegga function of the limma package 14 .  www.nature.com/scientificdata www.nature.com/scientificdata/ Data visualization. Ternary plot visualization was performed as previously described 15 . Ternary plots position cells according to the proportion of basal, LP-or ML-positive signature genes expressed by that cell and were generated using the vcd package 16 . The t-SNE visualization for all the integration analyses were generated using the RunTSNE function in Seurat with a random seed of 2018 for reproducibility. Diffusion plots were generated using the destiny package 17 . Multi-dimensional scaling (MDS) plots were created with edgeR's plotMDS function. Log2-CPM values for each gene across cells were calculated using edgeR's cpm function with a prior count of 1. Heat maps were generated using the pheatmap package. Log2-CPM values were standardized to have mean 0 and standard deviation 1 for each gene before producing the heat maps, after which genes and cells were clustered by the Ward's minimum variance method 18 . Bulk RNA-seq data and differential expression analysis. RNA-seq experiments were performed to obtain signature genes of basal, luminal progenitor (LP), mature luminal (ML) and stromal cell populations. Epithelial cells for basal, LP, and ML populations were sorted from eight independent patients and stroma from five patients. For one particular patient, samples were collected from both left and right breast for each of the four cell populations. For another patient, the ML cell population was collected twice. The complete RNA-seq data contains 9 basal, 9 LP, 10 ML and 6 stroma samples. RNA-seq libraries were prepared using Illumina's TruSeq protocol and were sequenced on an Illumina NextSeq 500.
Reads were aligned to the hg38 genome using Rsubread v1.5.3 19 . Gene counts were quantified by Entrez Gene IDs using featureCounts and Rsubread's built-in annotation 20 . Gene symbols were provided by NCBI gene annotation dated 29 September 2017. Immunoglobulin genes as well as obsolete Entrez Ids were discarded. Genes with count-per-million above 0.3 in at least 3 samples were kept in the analysis. TMM normalization was performed to account for the compositional biases between samples.
Differential expression analysis was performed using limma-voom 21 . Patients were treated as random effects and the intra-patient correlation was estimated by the duplicateCorrelation function in limma. Pairwise comparisons between the four cell populations were performed using TREAT with a fold change threshold of 1.5 22 . An FDR cut-off of 0.05 was applied for each comparison. Genes were considered as signature genes for a particular cell type if they were upregulated in that cell type in all the pairwise comparisons. The analysis yielded 515, 323, 765, and 1094 signature genes for basal, LP, ML, and stroma, respectively. In this submission we provide gene symbols of the signature genes as an R data file and R code to reproduce the bulk RNA-seq analysis 9 Differential abundance analysis. Differential abundance analyses were performed to examine the differences in cell cluster frequencies between pre-menopause and post-menopause groups in normal breast tissue www.nature.com/scientificdata www.nature.com/scientificdata/ micro-environment. Quasi-multinomial and quasi-binomial generalized linear models were used in order to account for the inter-patient variability. The numbers of cells under all the clusters from each individual donor were counted and used as the response variable in the model. The glm function of the stats package was used to fit the cell numbers against cell clusters, donors, plus a cluster-menopausal interaction term. The quasi-Poisson family was used in the glm function.
A quasi-multinomial F-test was performed to test for differences in cluster frequencies across all the clusters between pre-and post-menopausal samples, which yielded a p-value of 0.007. To test for cluster frequency differences for each individual cluster, we compared the cell numbers of that cluster with the aggregated cell numbers of all the other clusters across all the donors. Quasi-binomial generalized linear models were fitted and quasi-binomial F-tests were performed for each cluster separately. The p-values are 0.040 and 0.032 for cluster 1 and cluster 2, respectively, indicating these two clusters have significantly different sizes between pre-and post-menopause conditions after accounting for inter-patient variability. Sizes are not significantly different for other clusters. The R code to reproduce the differential abundance analysis is provided in the files NormEpi.R and NormTotal.R (Table 2). copy number variation analysis. Copy number variation (CNV) analysis was performed using inferCNV of the Trinity CTAT Project (https://github.com/broadinstitute/inferCNV), which compares gene expression intensity across genomic locations in the tumor or lymph-node samples with those in a normal reference sample. The single-cell RNA expression profile of a normal breast total cells sample (N-0372-total) was adopted as a reference for all the CNV analyses presented in the ScBrAtlas study. The results of each CNV analysis were visualized in a heatmap, which showed the relative expression intensities of the tumor samples with respect to the normal reference. For ease of visualization, cells from the same patient within the same cluster were grouped into a single column block, and only the blocks containing more than 100 cells were used in the heatmap. All the column blocks were assigned an equal width in each of the heatmap. The column block annotation of all the CNV heatmaps in this study is available as part of the Figshare deposition, indicating which clusters in which samples were classified as normal or tumor 9 .

Data Records
Cell Ranger genewise read counts for the 69 scRNA-seq profiles, prior to quality filtering, are available as GEO series GSE161529 7 . Quality filtering thresholds, downstream R data objects storing cell cluster identities and associated R code are available from Figshare 9 . Specific files available from Figshare are listed in Table 2. www.nature.com/scientificdata www.nature.com/scientificdata/ The bulk RNA-seq genewise read counts are available as GEO series GSE161892 23 . The cell-type signature genes generated from the bulk RNA-seq and associated R code are available from Figshare 9 .

Technical Validation
Technical quality of the 10X single-cell transcriptomic datasets was assessed by examining the number of mapped reads and the number of detected genes (genes with at least one read count mapped to it) for all cells across all the samples (Fig. 2a,b).
Quality control was performed to remove cells of low quality. Cells with a high proportion of mitochondrial reads or a low number of detected genes were removed. For each sample, an upper limit of library size was also used in combination with an upper limit of number of detected genes to remove potential multiplets. The proportion of cells retained after filtering is 82.2% across all 69 samples, indicating good data quality (Fig. 2c).
Technical quality of the bulk RNA-seq data was assessed using MDS and biological coefficient of variation (BCV) plots (Fig. 3).

Usage Notes
The code provided may be run using the free R programming environment with Bioconductor and Seurat R software packages https://www.r-project.org. The RDS files may be read using R's readRDS() function. The Seurat objects allow readers to use and extend the results of the major analyses conducted as part of the ScBrAtlas study. Cell barcodes and Seurat cell clustering information are stored in the meta.data component of each Seurat object.  www.nature.com/scientificdata www.nature.com/scientificdata/ code availability The R code files provided on Figshare contain complete the analyses of the ScBrAtlas study 9 (Table 2). Code files are also available from the GitHub repository https://github.com/yunshun/HumanBreast10X. All the bioinformatics analyses were performed in R 3.6.1 on x86_64-pc-linux-gnu (64-bit) platform, running under CentOS Linux 7. The following software packages were used for the analyses: Seurat v3. 1